In this example we’re going to use the built in faithful dataset which contains two columns: eruptions and waiting. We want to determine of the waiting times between Old Faithful eruptions are autocorrelated.

\[ H_0: \text{there is no correlation} \\ H_a: \text{there is a correlation} \]

First, let’s look at the data for Old Faithful waiting times.

waiting <- faithful$waiting
plot(waiting, type = "b")

Now, to estimate autocorrelation, we can do the following:

(observed_cor <- cor(waiting[-length(waiting)], waiting[-1]))
[1] -0.5397898

Now, in order to determine if this value is significant, we need some way to determine a p value. Once again, we can rely on permutation testing to build a distribution of values under the null hypothesis and then compare our observed value to that distribution.

n_permutations <- 10000
results <- replicate(n_permutations, {
  new_waiting <- sample(waiting)
  cor(new_waiting[-length(new_waiting)], new_waiting[-1])
})

Now that we have results, which contains 10000 autocorrelation values under the null hypothesis, we can examine its distribution.

plot(density(results))

Now, let’s look at our observed_cor value in light of these results.

plot(density(results), xlim = c(observed_cor, max(results)))
abline(v = observed_cor, col = "red")

Given our observed_cor and the values contained in results, we can caluclate a p value and associated confidence interval as follows:

(p_value <- mean(abs(results) >= abs(observed_cor)))
[1] 0
ci <- p_value + c(-1, 1) * qnorm(.975) * sqrt(p_value * (1 - p_value) / n_permutations)
c(lower = ci[1],
  p_value = p_value,
  upper = ci[2])
  lower p_value   upper 
      0       0       0
LS0tCnRpdGxlOiAiT2xkIEZhaXRoZnVsIC0gVGltZSBTZXJpZXMgQXV0b2NvcnJlbGF0aW9uIgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCgpgYGB7ciBzZXR1cCwgaW5jbHVkZSA9IEZBTFNFfQprbml0cjo6b3B0c19jaHVuayRzZXQoY2FjaGUgPSBUUlVFKQpgYGAKCkluIHRoaXMgZXhhbXBsZSB3ZSdyZSBnb2luZyB0byB1c2UgdGhlIGJ1aWx0IGluIGBmYWl0aGZ1bGAgZGF0YXNldCB3aGljaCBjb250YWlucwp0d28gY29sdW1uczogYGVydXB0aW9uc2AgYW5kIGB3YWl0aW5nYC4gV2Ugd2FudCB0byBkZXRlcm1pbmUgb2YgdGhlIHdhaXRpbmcgdGltZXMKYmV0d2VlbiBPbGQgRmFpdGhmdWwgZXJ1cHRpb25zIGFyZSBhdXRvY29ycmVsYXRlZC4KCiQkCkhfMDogXHRleHR7dGhlcmUgaXMgbm8gY29ycmVsYXRpb259IFxcCkhfYTogXHRleHR7dGhlcmUgaXMgYSBjb3JyZWxhdGlvbn0KJCQKCkZpcnN0LCBsZXQncyBsb29rIGF0IHRoZSBkYXRhIGZvciBPbGQgRmFpdGhmdWwgd2FpdGluZyB0aW1lcy4KCmBgYHtyfQp3YWl0aW5nIDwtIGZhaXRoZnVsJHdhaXRpbmcKcGxvdCh3YWl0aW5nLCB0eXBlID0gImIiKQpgYGAKCk5vdywgdG8gZXN0aW1hdGUgYXV0b2NvcnJlbGF0aW9uLCB3ZSBjYW4gZG8gdGhlIGZvbGxvd2luZzoKYGBge3J9CihvYnNlcnZlZF9jb3IgPC0gY29yKHdhaXRpbmdbLWxlbmd0aCh3YWl0aW5nKV0sIHdhaXRpbmdbLTFdKSkKYGBgCgpOb3csIGluIG9yZGVyIHRvIGRldGVybWluZSBpZiB0aGlzIHZhbHVlIGlzIHNpZ25pZmljYW50LCB3ZSBuZWVkIHNvbWUgd2F5IHRvIApkZXRlcm1pbmUgYSBwIHZhbHVlLiBPbmNlIGFnYWluLCB3ZSBjYW4gcmVseSBvbiBwZXJtdXRhdGlvbiB0ZXN0aW5nIHRvIGJ1aWxkIGEgCmRpc3RyaWJ1dGlvbiBvZiB2YWx1ZXMgdW5kZXIgdGhlIG51bGwgaHlwb3RoZXNpcyBhbmQgdGhlbiBjb21wYXJlIG91ciBvYnNlcnZlZAp2YWx1ZSB0byB0aGF0IGRpc3RyaWJ1dGlvbi4KCmBgYHtyfQpuX3Blcm11dGF0aW9ucyA8LSAxMDAwMApyZXN1bHRzIDwtIHJlcGxpY2F0ZShuX3Blcm11dGF0aW9ucywgewogIG5ld193YWl0aW5nIDwtIHNhbXBsZSh3YWl0aW5nKQogIGNvcihuZXdfd2FpdGluZ1stbGVuZ3RoKG5ld193YWl0aW5nKV0sIG5ld193YWl0aW5nWy0xXSkKfSkKYGBgCgpOb3cgdGhhdCB3ZSBoYXZlIGByZXN1bHRzYCwgd2hpY2ggY29udGFpbnMgMTAwMDAgYXV0b2NvcnJlbGF0aW9uIHZhbHVlcyB1bmRlciB0aGUKbnVsbCBoeXBvdGhlc2lzLCB3ZSBjYW4gZXhhbWluZSBpdHMgZGlzdHJpYnV0aW9uLgoKYGBge3J9CnBsb3QoZGVuc2l0eShyZXN1bHRzKSkKYGBgCgpOb3csIGxldCdzIGxvb2sgYXQgb3VyIGBvYnNlcnZlZF9jb3JgIHZhbHVlIGluIGxpZ2h0IG9mIHRoZXNlIHJlc3VsdHMuCgpgYGB7cn0KcGxvdChkZW5zaXR5KHJlc3VsdHMpLCB4bGltID0gYyhvYnNlcnZlZF9jb3IsIG1heChyZXN1bHRzKSkpCmFibGluZSh2ID0gb2JzZXJ2ZWRfY29yLCBjb2wgPSAicmVkIikKYGBgCgpHaXZlbiBvdXIgYG9ic2VydmVkX2NvcmAgYW5kIHRoZSB2YWx1ZXMgY29udGFpbmVkIGluIGByZXN1bHRzYCwgd2UgY2FuIGNhbHVjbGF0ZQphIHAgdmFsdWUgYW5kIGFzc29jaWF0ZWQgY29uZmlkZW5jZSBpbnRlcnZhbCBhcyBmb2xsb3dzOgpgYGB7cn0KKHBfdmFsdWUgPC0gbWVhbihhYnMocmVzdWx0cykgPj0gYWJzKG9ic2VydmVkX2NvcikpKQpgYGAKCmBgYHtyfQpjaSA8LSBwX3ZhbHVlICsgYygtMSwgMSkgKiBxbm9ybSguOTc1KSAqIHNxcnQocF92YWx1ZSAqICgxIC0gcF92YWx1ZSkgLyBuX3Blcm11dGF0aW9ucykKYyhsb3dlciA9IGNpWzFdLAogIHBfdmFsdWUgPSBwX3ZhbHVlLAogIHVwcGVyID0gY2lbMl0pCmBgYAoK